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Abstract 

The evolution of negative streamers during electric breakdown of a non-attaching 
gas can be described by a two-fluid model for electrons and positive ions. It con- 
sists of continuity equations for the charged particles including drift, diffusion and 
reaction in the local electric field, coupled to the Poisson equation for the electric 
potential. The model generates field enhancement and steep propagating ionization 
fronts at the tip of growing ionized filaments. An adaptive grid refinement method 
for the simulation of these structures is presented. It uses finite volume spatial dis- 
cretizations and explicit time stepping, which allows the decoupling of the grids for 
the continuity equations from those for the Poisson equation. Standard refinement 
methods in which the refinement criterion is based on local error monitors fail due 
to the pulled character of the streamer front that propagates into a linearly unstable 
state. We present a refinement method which deals with all these features. Tests 
on one-dimensional streamer fronts as well as on three-dimensional streamers with 
cylindrical symmetry (hence effectively 2D for numerical purposes) are carried out 
successfully. Results on fine grids are presented, they show that such an adaptive 
grid method is needed to capture the streamer characteristics well. This refinement 
strategy enables us to adequately compute negative streamers in pure gases in the 
parameter regime where a physical instability appears: branching streamers. 
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1 Introduction 



When non-ionized or lowly ionized matter is exposed to high electric fields, non-equilibrium ionization 
processes, so-called discharges, occur. They may appear in various forms depending on the spatio- 
temporal characteristics of the electric field and on the pressure and volume of the medium. One 
distinguishes the dark, glow or arc discharges that are stationary, and transient non-stationary phe- 
nomena such as leaders and streamers. We will focus here on streamers, that are growing filaments of 
plasma whose dynamics are controlled by highly localized and nonlinear space charge regions. 

Streamers occur in nature as well as in many technical applications. They play a role in creating the 
path of sparks and lightning [1] and they are believed to be directly observable as the multiple channels 
in so-called sprite discharges. These huge, lightning related discharges above thunderclouds attract 
large research effort since their first observation in 1989 [2,3,4,5]. Because of the reactive radicals they 
emit, streamers are used for the treatment of contaminated media like exhaust gasses [6,7], polluted 
water [8,9] or biogas [10]. More recently, efforts have been undertaken to improve the flow around 
aircraft wings by coupling space charge regions to gas convection [11]. 

Streamers can be either cathode directed or anode directed; the charge in their head is then positive 
or negative, respectively. Positive streamers propagate against the drift velocity of electrons, therefore 
they need additional (and not well known) mechanisms like nonlocal photoionization or background 
ionization. Negative streamers in simple non-attaching gases like nitrogen or argon, on the other hand, 
can be described by a minimal model with rather well-accessible parameters [12]. We therefore focus 
on this case. 

The minimal streamer model for negative streamers is a continuum approximation for the densities of 

the electrons and positive ions with a local field-dependent impact ionization term and with particle 
diffusion and particle drift in the local electrical field. As space charges of the streamer change the 
field, and the field determines drift and reaction rates of the particles, the model is nonlinear, and 
steep ionization fronts between ionized and non-ionized regions emerge dynamically. 

After the first claim [13] that a streamer filament within this deterministic continuum model in a 
sufficiently high field can evolve into an unstable state where it branches spontaneously, streamer 
branching was observed in more simulations [14,15,16]. While the physical nature of the Laplacian 
instability was elaborated in simplified analytical models [13,17,18,19,20], other authors challenged the 
accuracy of the numerical results [21,22,23,24]. Based on purely numerical evidence, their questions 
were justified, as all simulations within the minimal model up to now were carried out on uniform 
grids, and the numerical convergence on finer grids could hardly be tested. 

Therefore, in the present paper, we present a grid refinement strategy for the minimal streamer model 
and discuss its results. This procedure allows us to test the numerical convergence of branching, and 
also to calculate streamers efficiently in longer systems without introducing too many grid points. 
Moreover, the simulations in [13,14] were performed on the same uniform grid for both the particle 
densities and the electric potential. As the Poisson equation for the electric potential has to be solved 
on the complete large non-ionized outer space, this grid had to be much larger than the actual domain 
in which the particles evolve. Furthermore, the streamer has an inner layered structure with very steep 
ionization fronts and thin space charge layers. An efficient refinement is therefore badly needed. 
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There is an additional complication: standard grid refinement procedures in regions with steep gradi- 
ents fail due to the pulled character of the streamer front. "Pulling" means that the dynamics and, in 
particular, the front velocity, is determined in the linearly unstable high field region ahead of the front 
rather than by the regions with the steepest density gradients [12,25]. The high field region where 
any single electron immediately will create an ionization avalanche, is generated by the approaching 
curved ionization front itself. 

In the paper, a refinement strategy dealing efficiently with all these specific problems is developed. It 
is based on physical insight on the one hand, by which the relevant region for the particle densities 
can be restricted, and on the knowledge of the importance of the leading edge on the other hand. The 
drift-diffusion equations for the particle densities and the Poisson equation for the electric potential 
are treated separately, allowing the solutions for particle densities and electric potential to adapt to the 
specific difficulties. This refinement procedure will be applied to problems in two spatial dimensions, 
or in three dimensions with cylindrical symmetry. 

The outline of this paper is as follows. First, in Section 2, we give a brief description of the model. 
In Section 3 the numerical discretizations are introduced and motivated. Section 4 discusses the 
refinement procedure. Section 5 contains a one-dimensional example in which some of the essential 
numerical difficulties with local grid refinements of pulled fronts are illustrated and discussed. Section 6 
deals with the performance of our refinement algorithm, and the results are presented in the Section 7. 
The final section contains a discussion of the results and conclusions. 



2 Hydro dynamic approximation for the ionized channel 



2.1 Drift, diffusion and ionization in a gas 



The essential properties of anode directed streamer propagation in a non-attaching gas like or Ar, 
can be analyzed by a fluid model for two species of charged particles (electrons and positive ions). It 
consists of continuity equations for the electron and positive ion number densities, Wg and n+. 



dn 

-V ■ (ne/XeE + DeVne) = Si , (2.1) 

^ (2.2) 

The electrons drift with a velocity and diffuse with a diffusion coefficient Df,. Here E is the local 
electric field and He the electron mobility. The ions can be considered to be immobile on the short 

time scales considered here, since their mobility is two orders of magnitude smaller than that of the 
electrons [26]. We remark that extending our algorithm to moving ions and eventually other species 
is rather straightforward, and not including this makes no difference for the algorithm itself. 

The source term Si represents creation of electrons and positive ions through impact ionization, and is 
the same for both species since charge is conserved diiring an ionization event. The impact ionization 
can be described with Townsend's approximation [27] (but any other local field dependence can be 
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inserted as well) 



Si = ne/Xe|E|a(|E|) = neMe|E|ao e-^o/l^l , (2.3) 

where ao is the ionization coefficient and Eq the threshold field for ionization. We do not include 
photoionization or recombination since, in the particular case of N2 and for the short time scales 
considered, these processes are negligible compared to impact ionization [22]. However, an ionization 
mechanism such as photoionization is essential for the development of positive streamers, which are 
therefore excluded from the present study. 

The electric field E is determined through Poisson's equation for the electric potential V, 

VV = -(ne-n+), E = -Vy, (2.4) 
where eo is the permittivity of free space, and e the electron charge. 



2.2 Dimensional analysis 



This model has been implemented in dimensionless form. The characteristic length and field scales 
emerge directly from Townsend's ionization formula (2.3) as Iq = Oq ^ and Eq, respectively. The 
characteristic velocity is then given as vq = UsEq, which leads to a characteristic timescale to = 
^o/^^o = lo/iPeEo). The characteristic diffusion coefficient then becomes Dq = I^/Iq. The number 
density scale emerges from the Poisson equation (2.4), uq = eoEo/elo. We use the values from [26,28] 
for f^g, Eq and ao in N2 at 300 K, which depend on the neutral gas density, 

380 cm2 4332 1 ^ 2 • 10^ V , , 

"0 - TiTTT^ — , ^0 - TiTTTirTT — • (2-5) 



(N/Nq) Vs' "~(7V/7Vo) cm' ^ ~ (iV/iVo) cm' 

Inserting these values in the characteristic scales we obtain, for molecular nitrogen at normal condi- 
tions, 

2.3-10-4 3-10-12 1 1.8 -10^ cm2 ^ ^ 

Here Nq and To are the neutral gas density and temperature under normal conditions. The dimen- 
sionless quantities are then defined as follows, 

k to no no Eq (Eolo) Dq 

For the diffusion coefficient we use the value given in [29], 0^=1800 cm^ s"^, which gives a dimen- 
sionless diffusion coefficient of 0.1 under normal conditions. 
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Inserting these dimensionless quantities into the continuity equations (2.1)-(2.4), we obtain 



= V • (ct£ + DVa) + a\£\ exp(-l/|5|) , (2.8) 

^ = a\£\expi-l/\£\), (2.9) 

V2<^ = (T-p, f = -V(/), (2.10) 

where a and p denote the dimensionless electron and ion number densities, respectively, r the di- 
mensionless time, £ the dimensionless electric field, (f) the dimensionless electric potential and D the 
dimensionless diffusion coefficient. 

We refer to the equations (2.8)-(2.10) as the minimal streamer model, since it contains all the basic 
physics needed for negative streamers in a non-attaching gas. 



2.3 Geometry and boundary conditions 



In narrow geometries, streamers frequently are growing from pointed electrodes, that create strong 
local fields in their neighborhood and a pronounced asymmetry between the initiation of positive and 
negative streamers [30]. On the other hand, in many natural discharges and, in particular, for sprites 
above thunderclouds [15], it is appropriate to assume that the electric field is homogeneous. Of course, 
dust particles or other nucleation centers can play an additional role in discharge generation, but in 
this paper we will focus on the effect of a homogeneous field on a homogeneous gas. 

The computational domain is limited by two planar electrodes. The model is implemented in a cylin- 
drical symmetric coordinate system {r,z) £ (0, L^) x (0, L^), such that the electrodes are placed 
perpendicular to the axis of symmetry {r=0), the cathode at 2; = and the anode aX z = L^. The 
boundary conditions for the electric potential read 

<P{r,0,T) = 0, </.(r,L„r) = <^o>0, |^(L„z,r)=0. (2.11) 
The background electric field then becomes 



fb = -|£b|ez = -7^e,, (2.12) 

where is the unit vector in the z— direction. The radial boundary at Lr is virtual, and only present 
to create a finite computational domain. In order for the boundary condition no to affect the solution 
near the axis of symmetry along which the streamer propagates, we need to place this boundary 
relatively far from the axis of symmetry. 

Throughout this article, we use a Gaussian initial ionization seed, placed on the axis of symmetry, at 
a distance z = Zf, from the cathode. 
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a(r, z, 0) = p(r, z, 0) = ao exp {L^S^^^ . (2.13) 

The maximal density (Jq of this seed, the radius i?b at which the density drops to 1/e of its maximal 
value, and the value of differ from case to case, and will be specified where needed. Furthermore, 
the use of a dense seed, in particular in low fields, accelerates the emergence of a streamer. We remark 
that the initial seed is charge neutral. 

The continuity equation for the electron density is second order in space, and therefore requires two 
boundary conditions for each direction in space. At r = and z = we use Neumann boundary 
conditions, so that electrons that arrive at those boundaries may flow out of or into the system, but in 
all cases discussed in this paper the streamer is too far from the boundary for this to happen. At the 
cathode, we impose either homogeneous Neumann or Dirichlet boundary conditions. In the first case 
we again allow for a net flux of particles through the boundary. Dirichlet conditions will only be used 
for a one-dimensional test in this paper. To recapitulate, the boundary conditions for the electrons 
read 



^(r,0,r) = 0, or a(r,0,T) =0, ^(r,L„r) = 0, -£{Lr, z,t) = . (2.14) 

We notice that, if Zb » Rb, the ionization seed is detached from the cathode, and this results in 
practice in a zero inflow of electrons. On the other hand, placing the seed near the cathode will result 
in an inflow of electrons. Varying the value of Zb will therefore enable us to investigate the influence 
of the inflow of electrons on the streamer propagation. 



3 Numerical discretizations 



In our numerical simulations we shall mainly consider the streamer model with radial symmetry, 
making it effectively two dimensional. To illustrate some of the difficulties and their solutions we will 
also deal with the one-dimensional case. 

In the cylindrically symmetric coordinate system introduced in the previous section, the equations (2.8)- 
(2.10) read 



da _1 d{raSr) _^ d{a£z) d _|_ jj'^^'^ ^ 151 I 

dr r dr dz r dr dr dz'^ ' 



(3.1) 



^ = a\£\e-y\^\, (3.2) 

2 , 1 d d(j) d'^<p 

^=rd-r^'d^^^d^ = ^-P- ^'-'^ 
The electric field £ = {£r,Sz)^ can be computed from the electric potential as 
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(3.4) 



The boundary conditions for this system have been treated in Sect. 2.3. 



3. 1 Spatial discretization of the continuity equations 

The equations will be solved on a sequence of (locally) uniform grids with cells 

= [{i - l)Ar,iAr] x [(j - l)Az, jAz], i = l...M^, j = 1...M, , 

where M,. = Lr/Ar and = Lz/Az arc the number of grid points in the r- respectively z-direction, 
and cell centers (rj, zj) = {{i — ^)Ar, {j — i)Az). We denote by aij and pij the density approximations 
in the cell centers. These can also be viewed as cell averages. The electric potential cpij and field 
strength \£\i,j are taken in the cell centers, whereas the electric field components are taken on the 
cell vertices. For the moment it is supposed that the electric field is known, its computation will be 
discussed later on. 

The equations for the particle densities are discretized with finite volume methods, based on mass 
balances for all cells. Rewriting the continuity equations (3.1)-(3.2) will result in the semi-discrete 
system 

(3.5) 

dr ~ '''' 

in which F = F"- + F'^. F"' and F'^ are the advective and diffusive electron fluxes through the cell 
boundaries, and Sij is the source term in the grid cell Cij. 

The discretization of the advective terms requires care. A first order upwind scheme as used in [22,24] 
appears to be much too diffusive [31], leading to a totally different asymptotic behavior on realistic 
grids. Moreover, it is expected that the numerical diffusion might over-stabilize the numerical scheme, 
thereby suppressing interesting features of the solutions. This explains why streamer branching is not 
seen by the authors of [22,24]. On the other hand, higher order linear discretizations lead to numerical 
oscillations and negative values for the electron density, that will grow in time due to the reaction term. 
This holds in particular for central discretizations [32]. The choice was made to use an upwind-biased 
scheme with flux limiting. This gives mass conservation and monotone solutions without introducing 
too much numerical diffusion. For the limiter we will take the Koren limiter function, which is slightly 
more accurate than standard choices such as the van Leer limiter function [32]. 

Denoting = max(u, 0) and v~ = min(f , 0) to distinguish upwind directions for the components of 
the drift velocity Vr = —£r S'ld Vz = Sz, the advective fluxes in the r- direction are computed by 
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in which 



V'W = max(0,inin(l,i + ^,e)) . (3.7) 



The advective fluxes in the vertical direction arc computed in the same way; the fact that the r- 
direction is radial is already taken care of in (3.5). Note that in regions where the solution is smooth 
we will have values of 9ij close to 1, and then the scheme simply reduces to the third-order upwind- 
biased discretization corresponding to tp^O) = ^ + ^6. In non-smooth regions where monotonicity is 
important the scheme can switch to first-order upwind, which corresponds to 'tp{6) = 0. 

The diffusive fluxes are calculated with standard second-order central differences as 

Finally, the reaction term Sij in (3.5) is computed in the cell centers as 

Sij = C7i,,|5|i,e-Vl£l«> . (3.9) 

Boundary values will be either homogeneous Dirichlet or homogeneous Neumann type. For example, 
for Dirichlet boundary conditions o" = for z = 0, we introduce virtual values [32] 

o'i.o = , cTj-i = -o-j,2 • (3.10) 
For Neumann boundary conditions dzcr = for 2; = we set 

Ci.O = Ci,! , CTi^-l = cri,2 ■ (3-11) 
These formulas follow from the approximations 

o-z{ri,zi) = —{a{ri,zi)-a{ri,zo)) + 0{Az ) 
3.2 Spatial discretization of the Poisson equation 



The electric potential (p is computed through a second-order central approximation of Eq. (2.10), and 
is defined at the cell centers: 
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_ _ (l>i+i,j - '2(t>i,j + 4>i-i,j , (t>i+i,j - 4>i-i,j , <l>i,j+i - + <t>i,j-i 



The electric field components are then computed by using a second-order central discretization of 
£ = — V(/>, they axe defined in the cell boundaries, 



The electric field strength is computed at the cell centers, therefore the components are first deter- 
mined in the cell centers by averaging the cell boundary values, and the electric field strength then 
becomes: 



l\/(^r,-y + £r,^,j)' + (^.,,-1 + ■ (3-15) 

We notice here that, discretizing V • £ with a second-order central scheme gives 

Therefore, the total current conservation, 

V ■ (^^ + a£ + DVa^ =0, (3.16) 
also holds on the level of the discretizations. 

3.3 Temporal discretization 



After the spatial discretization, the system of equations (3.5) can be written in vector form as a system 
of ordinary differential equations, 

dcr ^, „^ 
^ = G(a,£), 

; (3.17) 

where the components of G and S are given by the spatial discretizations in (3.5). The electric field 
£ and the field strength \£\ are computed from given cr, p by discretized versions of (3.3) and (3.4), 
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discussed in Sect. 4.4. Therefore the full set of semi-discrete equations actually forms a system of 
differential-algebraic equations. 

The particle densities are updated in time using the explicit trapezoidal rule, which is a two-stage 

method, with step size At. Starting at time = nAr from known particle distributions a'^{r, z) ~ 
a{r,z,T^), p'^{r,z) « p{r,z,T"') and known electric field E^{r,z) « £{r,z,T^), the predictors a"'^^ 
and p""*"^ for the electron and ion densities at time t'^+^ are first computed by 

-n+i ^ a" + ATG(a",£"), 

(3.18) 

-n+i ^ p" + Ar5(a",£"), 



After this the Poisson equation (3.3) is solved with source term 0""+^ — /f'^^, leading to the electric 
field at this intermediate stage by Eq. (3.4). The final densities at the new time level r"'^^ are 

then computed as 

a ^ = a + —U{a ,t ) + —G[a ^ ,t ), 

(3-19) 

after which the Poisson equation is solved once more, but now with the source term o"""'"^ — p""*"-^, giving 
the electric field 5"+^ induced by the final particle densities at time r""*"^. This time discretization is 
second-order accurate, which is in line with the accuracy of the spatial discretization. 

The use of an explicit time integration scheme implies that the grid spacing and time step should obey 
restrictions for stability. For the advection part we get a Courant-Friedrichs-Lewy (CFL) restriction 

At At 

maxf^-j hmaxf^^— — < i^a , (3.20) 

Ar Az 

and the diffusion part leads to 



^|S+^^<-<^- (3-21) 

Actually, to be more precise, a combination of (3.20), (3.21) should be considered. However, in our 
simulations, condition (3.20) will dominate and the time step will be chosen well inside this constraint. 

For the first order upwind advection scheme combined with a two-stage Runge-Kutta method, the 
maximal Courant number is [32] 1/1 = 1, while for the third order upwind scheme it is = 0.87. The 
second order central discretization demands a maximal Courant number = 0.5. 

A third restriction for the time step comes from the dielectric relaxation time. The ions are considered 
to be immobile, which leaves us with the following time step restriction in dimensional units [33], 
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At < , (3.22) 

efie maxne 

where we refer to the previous section for the meaning of these quantities. If we apply the dimensional 
analysis of the previous section, we obtain the time step restriction in dimensionless units, 



Ar < - . (3.23) 

In practice, it appears that this restriction is much weaker than that for the stability of the numerical 
scheme. 

The choice for an explicit time stepping scheme was made after tests with VLUGR [34], a local 
refinement code that uses -computationally much more intensive- implicit schemes (BDF2). These 
tests showed that these implicit schemes also needed relatively small time steps to obtain accurate 
solutions, so that in the end an explicit scheme would be more efficient in spite of stability restrictions 
for the time steps. Moreover the use of an explicit scheme allows us to decouple the computation of the 
particle densities from that of the electric potential and electric fields. With a fully implicit scheme all 
quantities are coupled, leading to much more complex computations and longer computation times. 

Another reason for preferring explicit time stepping is monotonicity. The solutions in our model 
have very steep gradients for which we use spatial discretizations with limiters to prevent spurious 
oscillations. Of course, such oscillations should also be prevented in the time integration. This has 
led to the development of schemes with TVD (total variation diminishing) or SSP (strong stability 
preserving) properties; see [35,32]. In contrast to stability in the von Neumann sense (i.e., L2-stability 
for linear (izcd) problems with (frozen) constant coefficients), such monotonicity requirements do not 
allow large step sizes with implicit methods of order larger than one. Among the explicit methods, 
the explicit trapezoidal rule is optimal with respect to such monotonicity requirements. 



3.4 Remarks on alternative discretizations 



The above combination of spatial and temporal discretizations provide a relatively simple scheme 
for the streamer simulations. The accuracy will be roughly (9(A,t^) + O(Ar^) in regions where the 
solution is smooth (also for the limited advcction discretization [32, p. 218]). In our tests, step size 
restrictions mainly originate from the advective parts in the continuity equations. The above scheme 
is stable and monotone (free of oscillations) for Courant numbers up to one, approximately. Usually 
we take smaller step sizes than imposed by this bound to reduce temporal errors. 

As mentioned before in Sect. 3.1 using a first-order upwind discretization for the advective term will 
usually give rise to too much diffusion, whereas second-order central advection discretizations lead to 
numerical oscillations and negative concentrations. 

Higher-order discretizations can certainly be viable alternatives. However, we then will have larger 
spatial stencils, which creates more difficulties with local grid refinements where numerical interfaces 
are created. The above discretization is robust and easy to implement. 

It is well known that limiting as in (3.6) gives some clipping of peak values in linear advection tests. 
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simply because the limiter does not distinguish genuine extrema from oscillations induced numerically. 
This can be avoided by adjusting the limiter near extrema, but in the streamer tests it was found 
that such adjustments were not necessary. In the streamer model the local extrema in each spatial 
direction are located in the streamer head, and the nonlinear character of the equations gives a natural 
steepening there which counteracts local numerical dissipation. 

In [28] a flux-corrected transport (FCT) scheme was used. The advantage of our semi-discrete approach 
is that is becomes easier to add source terms without having to change the simulation drastically. 
Also the transition from ID discretizations to 2D or 3D becomes straightforward conceptually; the 
implementations for higher dimensions are still difficult, of course, in particular for the Poisson equa- 
tion. Moreover, in [31] comparisons of the FCT scheme with a scheme using the van Leer limiter 
(which is closely related to the limiter used in our scheme) show that the FCT scheme, in contrast to 
the limited scheme, gives somewhat irregular results for simple advection tests in regions with small 
densities. In the leading edge the densities decay exponentially, and we do not want such irregularities 
to occur. 



4 The adaptive refinement procedure 

4-1 The limitations of the uniform grid approach 

Up to now, all simulations that have been carried out on the minimal streamer model were performed 
on a uniform grid [26,28,13,14]. However the use of a uniform grid on such a large system has several 
limitations. 

The first limitation is the size of the system. In [13,14] the simulations were performed in a radially 
symmetric geometry with 2000 x 2000 = 4 • 10® grid points. Since the number of variables is of at least 
10 per grid point (the electron and ion densities both at old and new time step, the electric potential, 
the electric field components and strength, and the terms containing the temporal derivatives of both 
electrons and ions), the total number of variables will be at least 40 • 10®. So, when computing in 
double precision (64 bits or 8-bytes values) the memory usage is in the order of several hundreds 
MB, depending on the compiler. Moreover, these simulations show that the streamer will eventually 
branch, and up to now there was no convergence of the branching time with respect to the mesh size. 
In order to investigate the branching, it would be necessary to rerun the simulation with a smaller 
grid size. Moreover it would be interesting to investigate larger systems. So from that point of view 
it is worth looking at an algorithm that would require much less memory usage. 

The second limitation conies from the Poisson equation, which has to be solved at each time step, 
and therefore requires a fast solver. For this we use the Fishpak routine, described in [36,37]. One 
of the major limitations of this routine is its ineptitude to deal accurately with very large grids, due 
to numerical instabilities with respect to round-off errors. Numerical experiments show that on an 
rrir X rriz grid with either or substantially larger than 2000, the FiSHPAK routine gives large 
errors due to numerical instabilities with respect to round-off errors. An illustration is presented in 
the appendix. 
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Fig. 1. An example of nested grids: ■ cell centers on with = 1, x cell centers on with l{2) = 2, 
+ cell centers on with Z(3) = 3, o cell centers on Q'^ with Z(4) = 3. 

It is necessary to develop some strategy to counteract these limitations. This will be done by choosing 
separate grids with suitable local refinements for the continuity equations (2.8)-(2.9) and the Poisson 
equation (2.10). 

4-2 General structure of the locally uniform refined grids 

Both the continuity equations (3.1) and the Poisson equation (3.3) are computed on a set of uniform, 
radially symmetric grids, that are refined locally where needed. Solving these equations separately 
rather than simultaneously allows the use of different sets of grids for each equation, thereby making 

it possible to decouple the grids for the continuity equation from those of the Poisson equation; grids 
can then be tailored for the particular task at hand. We emphasize that it is the use of an explicit 
time stepping method that allows us to decouple the grids. 

In both cases the equations are first solved on a coarse grid. This grid is then refined in those regions 
where a refinement criterion - which will be treated in more detail below - is met. These finer grids can 
be further refined, wherever the criterion is still satisfied. Both the grids and the refinement criteria 
may be different for the continuity equations on the one hand and the Poisson equation on the other 
hand, but the general structure of the grids is the same for both type of equations. It consists of a 
series of nested grids Q,^, k being the grid number, with level l{k). This level function gives the mesh 
width of a grid, ^(1) = 1 corresponding to the coarsest grid fi^ with mesh width Ar'^ and Az"^. A 
certain grid will have a mesh twice as fine as the grid one level coarser - which we will call its parent 
grid - so that the mesh widths of a grid with level I become Ar' = Ar'^/2^~^ and Az^ = Az'^/2''~^. 

Fig. 1 shows an example of four nested grids on three levels. We will denote a quantity u on a certain 
grid ri*^ with level I = l{k) as uf^ = u{rl,Zj). All grids are characterized by the coordinates of 
their corners relative to the origin of the system (on the axis of symmetry at the cathode), so that 
rl = {i- 1/2) Ar', = {j - 1/2) Az^. 

Throughout this article, the grids for the Poisson equation are denoted by Q, those for the continuity 
equations by TC. For the continuity equations, the grids are structured as follows: we determine all 
grid cells of a grid Ti^ at a certain level I that, following some refinement monitor (we will give more 
details later), have to be refined. This results in one or more sets of connected finer grid cells. As a 
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first step, we chose the finer child grids {Ti'"''^} of the parent gridH} to cover the smallest rectangular 
domain enclosing each of these sets of connected grid cells. Although this first rudimentary approach 
gave some gain in computational time, it lead to a large number of unnecessary fine grid cells, due to 
the curved nature of the ionization front. So instead, we divide such unnecessary large rectangles into 
smaller rectangular patches, thereby limiting the number of fine grid cells. Obviously, the union of all 
these child grids {H}^^} should contain all the grid cells indicated by the monitor. For programming 
reasons, the rectangular child grids {Ti}^^} were chosen to all have an equal, arbitrary number Mq of 
grid points in the radial direction. Using a large value for Mq leads to unnecessary large child grids, 
using a small value would make no sense because each grid requires the storage of boundary values. 
Since the rectangular grid structure is preserved, this could be relatively easily be implemented in the 
existing code, and using a suitable value of Mq lead to a significant gain in computational time and 
memory. 

To compute the electric potential, however, such a grid distribution is not appropriate. This comes 
from our use of a fast Poisson solver, which computes, on a rectangular grid, the potential induced 
by a space charge distribution on that same rectangular grid. The solution of the Poisson equation 
is not, however, determined locally, and these non local effects arc not accounted for properly if we 
compute the potential on smaller rectangles like the continuity equation. So in the case of the Poisson 
equation we use the same grid structure as we first did for the continuity equation: we determine, 
using some refinement monitor, all grid cells of a grid at a certain level / that have to be refined, 
and the finer child grid {^'+^} of the parent grid is set to cover the smallest rectangular domain 
enclosing all these grid cells. 

In what follows, to make the distinction between the indices on a coarse and on a fine grid, we use 
capital indices I and J for the coarse grid and small indices i and j for the fine grid. Notice that 
a coarse grid cell with the cell center {rj,zj) contains four finer cells with centers {ri,Zj), {ri,Zj^i), 
{ri+i,Zj), (rj+i, Zj+i), with i = 2/ - 1 and j = 2J - 1. 

Let us now discuss in more details the refinement algorithm for the continuity equations and that for 
the Poisson equation, and the benefits of decoupling the grids for both equations. 



4-3 Refinement scheme for the continuity equation 

Let us assume that the particle distributions and electric field at time r" are known on a set 5" of m 
rectangular grids H"''^ with level I = l{k), 1 < k < m, as shown in Fig. 2. 

Then, using the explicit time stepping method introduced in Sect. 3.3, the particle distributions at 
time r"+^ can be computed on all the grids of S"' (Fig. 2b) . Now the new set /S""*"^ of nested grids 
that is best suitable for the solution at r"+^ has to be found, as in Fig. 2c. 

Moreover, the computational domain for the continuity equations can be reduced substantially by 
the following physical consideration: our model is a fluid model based in the continuum hypothesis, 
which is not valid anymore if the densities are below a certain threshold, that is taken as 1 mm~^. In 
nitrogen at atmospheric pressure, this corresponds roughly to a dimensionless density of 10^^^. The 
regions where all densities are below this threshold is ignored. Therefore, after each time step the 
densities below this threshold are set to zero. The computational domain for the continuity equations 
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for the next time step is then taken as the region where a or p are strictly positive. In view of our 
two-stage Runge-Kutta time stepping we use a four point extension of this domain in all directions. 



4.3.1 Restriction of fine grid values to a coarse grid 

At first, when a grid Ti"''^ at level L = L{K) contains a finer grid H"''^ at level Z = L + 1, the particle 
densities on 7^"'^ are restricted to W"'^ in such a way that the total mass of each particle species is 
conserved locally. Prom Fig. 1 it can be seen that one cell of H"'^ contains four cells on H"''^, and the 
mass conservation implies that for each particle species the total mass in the coarse grid cell is equal to 
the sum of the masses in the finer grid cells, which, taking into account the cylindrical geometry of the 
cells, translates in the following restriction formula = Res{U^) for the grid functions = {ufj} 
and = {<■}, 



21 2J 

Ufj = Res{U^)ij=-j; E ^-«-n, (4.1) 

m=2/-ln=2J-l 



in which u stands for either the electron or the ion density. This restriction step is carried out because 
time stepping on a too coarse grid may lead to erroneous values, which are now overwritten by the 
better restricted values. The r-factors account for the mass distribution in the radial cells in cylindrical 
symmetry. 




Fig. 2. (a): Contour line of the solution and the set of rectangular computational grids, both at r". (b): 
Contour line of the solution at t"+^ and computational grids at r"^. (c): Contour line and computational 
grids at r""^^, the dashed lines are the grids at r". 
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4-3.2 Refinement criterion: curvature monitor 

It is now possible to find the regions where the grids should be refined at The decision whether a 
finer grid should be used on a certain region is made with a relative curvature monitor. This monitor 
is defined on a grid with level I as the discretization of 



Mu{r\z^) = (Ar 



l\2 



19 du^ 
r dr dr' 



Z\2 



(4.2) 



Although this expression does not provide an accurate estimate of the discretization errors, it does 
give a good indication of the degree of spatial difficulty of the problem [38] . It is applied first to the 
electron density o", since that is the quantity which advects and diffuses, and therefore the quantity 
in which some discretization error will appear. The monitor is also applied to the total charge density 
a — p since this is the source term of the Poisson equation, and the accuracy of the solution of the 
Poisson equation is of course dependent on the accuracy of the source term. The monitor is taken 
relative to the maximum value of each quantity, and the refinement criterion through which the need 
to refine a certain grid Ti^ with level / then reads: 



refine all grid cells i,j where — ^ — " > e^, u = u, u — p (4.3) 

max 

''J 

in which and e^^_p are grid-dependent refinement tolerances that will further be discussed in Sect. 5 
and Sect. 6. 

Now, starting from the coarsest grid, the monitors are computed by approximating them with a second 
order central discretization, which determines the regions that should be refined. For this set of finer 
grids again the regions to be refined arc computed, and so on either until the finest discretization 
level is reached, or until the monitor is small enough on every grid. Now the new set of nested grids 
gn+i j^g^g ]-,ggj^ determined, but the particle densities are still only known on the set S"^ (in Fig. 2c 
this means that we have to convert the density distributions on the dashed grids to distributions on 
the solid rectangles). 

Criteria such as (4.2)-(4.3) are common for grid refinements. As we shall see in experiments, it will 
be necessary to extend the refined regions to include (a part of) the leading edge of a streamer. This 
leading edge is the high field region into which the streamer propagates, and where the densities decay 
exponentially. This modification, due to the pulled front character of the equations [25], is essential 
for the front dynamics to be well captured, and is a major new insight for the simulation of streamers, 
and more generally, of any leading edge dominated problem. It is discussed in more detail in Sect. 5. 



4-3.3 Grid mapping 

The refinement monitor gives a new grid distribution {H""'"^''^} at time r""*"^. In the following we 
consider mapping functions used to map the solution at r"^^ on the previous grid distribution {"H"'*^} 
on the new grid distribution {H""^^'^]. There are three possible relations between a grid T^""*"^'*^ at 
time r^'^^ and the older grids {W"''^}, as illustrated in Fig. 3. 



16 



Fig. 3. Three grids at time r" (solid lines) and a grid at next time step r"^^ (enclosed by the dashed 
lines) with the same level as W"'^. In the vertically striped region the new grid coincides with a finer grid 
at the previous time step, in the horizontally striped region only a coarser grid existed at r", and in the 
region that is not filled a grid at the same level existed at r". 

First it is possible that a fine grid H"'*^ at the previous time step now is covered by a coarser grid 
'j^n+i,K ^g^g ^jjg vertically striped region of Fig. 3). Then the values of the densities on the new grid 
are computed by the restriction (4.1): 



where U again stands for the set of grid values of either electron- or ion densities. 

Secondly, there is the possibility that (part of) the new grid already existed at previous time step, in 
which case there is no need for projecting the density distributions from one grid to the other. 

Finally, it may occur that the new grid lies on a region where only a coarser grid n""^ existed at 
the previous time level (as in the horizontally striped region in Fig. 3). Then we have to prolongate 
the coarse grid values to the new fine grid. One main consideration in the choice of the 

prolongation is the conservation of charge in the discretizations. For simplicity, we will first consider a 
one-dimensional prolongation, which is then easily extended to more dimensions. We consider coarse 
grid values = {uf}. Then a mass conserving interpolation for values = {uf} on a grid twice 



ljn+l,K ^ Ilf,s{ljn,k^^ ^n+l,K p ^n,k 



(4.4) 



as fine is, 



u. 



(4.5) 



which obviously implies mass conservation, Az^uf + Az^uf_^_-^ = Az^Uf. Using a three-points stencil, 
the coefficient Dj, such that the interpolation is second order accurate, can be written as 



Di = l{Ui^i + Ui+i). 



(4.6) 
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In the case of a three-dimensional geometry with axial symmetry, the above interpolation is applied 
on in the ^-direction, and on RU^ in the radial direction, which ensures mass conservation on 
such a system. 

Remark. This mass conserving interpolation might lead to new extrema or negative values. That could 
be prevented automatically by limiting the size of the slopes. However, in our simulations this turned 
out not to be necessary. 

Finally, we need boundary values for all grids. On the coarsest grid these simply follow from the 
discretization of the boundary conditions (2.11)-(2.14), as explained in Sect. 3.1. On finer grids they 
are interpolated bi-quadratically from coarse grid values. The interpolation error is then third order, 
and therefore smaller than the discretization error, which would not be the case with a bi-linear 
interpolation. The quadratic interpolation is derived using the Lagrange interpolation formula to find 
jjk _ ^yk.y from the coarse grid values = {ufj}. The interpolated value becomes 

Uij = Int{U^)ij = -j^ ^ ^ Tj^p uf+p,J+q n 1{K) _ 1{K) n 1{K) _ 1{K) ' ('^■'^) 

p=-lq=-l P^p I^I+p ^I+P Qj^q Zj+q ^J+Q 

P=-l Q=-l 

in which I = {i + i mod 2)/2 and J = {j + j mod 2)/2. We notice that the interpolated quantity 
again is not a but ra because of the cylindrical coordinate system. 

In our algorithm, mass conservation at the grid interfaces will simply be ensured by matching the 
fluxes of the fine and coarse grids. 

4.3.4 Flux conservation at grid interfaces 

The mapping from one grid to the other must take into account a flux correction on grid interfaces, 
in order to ensure mass conservation. This correction yields that the total mass going through one 
coarse grid cell must be the same as the sum of total mass going through two flne grid cells. Taking 
into account the cylindrical geometry of the system, the fluxes through coarse grid cells become, in 
terms of the fine grid fluxes, 

-^rf/+l/2,J = -^p:(K;i+l/2,j + ^r;i+l/2J+l) ' 
^z;I,J+l/2 = /^^L^j i4Fz;i,j+l/2 + 4+lF^;i+l,j+l/2) ' 

where i = 21 — 1 and j = 2J — 1. 

4.4 Refinement scheme for the Poisson equation 

The refinement strategy for the computation of the electric potential is also based on a static regridding 
approach, where the grids are adapted after each complete time step. In this case however, the 
refinement criterion is not based on a curvature monitor but on an iterative error estimate approach. 



18 







-o 



,1 



Fig. 4. The nine-points stencil used to map coarse grid values of the electric potential (f)f} onto the fine 
grid, thus obtaining cf>lj, 4>\+ij< ^'^'^ The cell centers of the coarse grid are marked with 

■ , those of the fine grids with o, and the boundary points of the fine grids are marked with x. 



The FiSHPAK routine will be used on a sequence of nested grids Q'^. The solution on a coarse grid will 
be used to provide boundary conditions for the grid on the finer level, which will in general extend 
over a smaller region. This approach is explained in full detail in [39]; here we will discuss the main 
features of the scheme. 

Starting from the (known) charge density distribution o" — p on a set of grids {H^}, the Poisson 
equation is first solved on the two coarsest grids and Q^, both covering the entire computational 
domain (0, L^) x (0, L^). The finest of these two grids is coarser or as coarse as the coarsest grid of 
{T^'^}. The densities should then first be mapped onto the coarse grids and Q"^, using the restriction 
formula (4.1). The source term of the Poisson equation is then known on these two coarse grids, and 
the Poisson equation is solved using a FiSHPAK routine that discretizes Eq. (3.3) using second-order 
differences: 



a, 



1,3 



+ 



2r|Ar' 



+ 



(4.9) 



in which I = l{m) is the level of grid Q"^. This system of linear equations is then solved using a cyclic 
reduction algorithm, see [40] . The details of FiSHPAK are described in [37] . The subroutine was used 
as a black box in our simulations. A comparison with iterative solvers, multigrid or conjugate gradient 
type, can be found in [41]. For the special problem we have here - Poisson equation on a rectangle 
- such iterative solvers are not only much slower than FiSHPAK, but they also require much more 
computational memory. 

As a next step, the coarse grid electric potentials 0^ on and 0^ on Q'^ are compared with each 
other, by mapping 0^ onto with a quadratic interpolation based on a nine-point stencil as shown 
in Fig. 4. 

For the prolongation that gives a continuous map of an electric potential 0^ = {(pfj} onto a finer 
grid Q'^ the following formula, based on a least square fit, is used: 
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(t)"'{r,z) = Pro{ (j)^ {r,z)) 

= (Pfj + cio(r - rf ) + CQi{z - zf) + cii(r - rf){z - zf )+ (4.10) 

+ C2o{r-T^f + c,2{z-z\-'f. 

For the values of the coefficients Cij we refer to [39]. This interpolation will primarily be used to 
generate boundary conditions for the finer level, as illustrated in Fig. 4 for the boundary points on 
the finer grid marked with x , but it will also be used for error estimation. The interpolation for the 
boundary conditions will be done such that there is a bias in the stencil toward the smooth region to 
enhance the accuracy. 

The error of the solution on a grid Q"^ is then estimated by the point-wise difference between the 
solution on that grid and the interpolation of the solution on the underlying coarser grid : 

~9T,i = - (^^« ^""kj ■ (4-11) 

The refinement criterion is as follows: 

refine all grid cells i,j where > (4-12) 

where is a small parameter that still has to be chosen. This leads to new and finer grids on which 
the whole procedure of mapping the charge densities onto it, solving the Poisson equation using the 
FiSHPAK routine and estimating the error is repeated, and so on either until the error is smaller than 
€0 everywhere or until the finest desired level is reached. Notice that, since the grids do not cover the 
whole computational domain anymore, the boundary conditions (2.11) do no longer hold. We now 
have Dirichlet boundary conditions on each boundary that lies in the interior of the computational 
domain, and they are computed by interpolating the electric potential on the finest grid that is known 
using Eq. (4.10). Using this third order interpolation formula implies that the error introduced by the 
interpolation is negligible compared to to discretization error, which is second order. With this method 
there is an upper bound for the maximum error on grid m with level I [39]: 

e"" <3maxg^ + {l-l)e^, (4.13) 

where g"^ = maxgfj as defined in Eq. (4.11). This means that the extra error due to the refinement 
can be made as small as wanted provided is taken small enough. Therefore, iteration, e.g. with 
defect corrections, is not needed. Inequality (4.13) is based on the assumption that the interpolation 
errors are negligible compared to the second-order discretization errors of the local problems, and 
therefore interpolant (4.10) was chosen to be of higher order. Although tacit smoothness assumptions 
are involved here, tests in [39] with strong local source terms did show that the errors are indeed well 
controlled by this nested procedure. 

We notice that this global error control is the reason for the choice of this particular refinement 
monitor for the Poisson equation. Such an error control does not hold for the continuity equations, 
for which we use the more suitable curvature monitor [34,38]. 

The electric field has to be known on the continuity grids 7i^, since it appears in the continuity 
equations (2.8)-(2.9). However it has to be computed from the electric potential that is only known 
on the Poisson grids, using Eq. (3.4). We consider the grid H'' with level l{k) on which the electric field 
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has to be known, and the finest potential grid G"^ with level l{m), which has a non-empty intersection 
with H''. 

There are three possible situations: the potential grid can be either coarser, as fine as, or finer than the 
continuity grid. If both grids have the same size {l{k) = l{m)), or if the continuity grid is coarser than 
the Poisson grid {l{k) < l{m)), then the field components are set directly with a second order central 
approximation of Eq. (3.4), using the neighboring points of the point where the field components need 
to be known. If the Poisson grid is coarser than the continuity grid {l{k) > l{m)), the electric field 
is first computed on the Poisson grid with a second order central approximation of Eq. (3.4), and 
then this field is interpolated to the continuity grid. The interpolation is performed with a piecewise 
bilinear approximation. 

4-5 Overall algorithm 

In previous sections, the refinement algorithms for the continuity equations and the Poisson equation 
have been treated separately. We here describe the overall algorithm. 

We start from known density distributions of the charged particles cr" and />" at a certain time r", 
on a set of grids H"''^. The electric field induced by the charges on the grids K"''^ is computed using 
the refinement method described in Sect. 4.4. 

The step size for the time integration is then set in such a way that the stability conditions (3.20)- 
(3.21) of both advection and diffusion discretizations are met on the finest grid. The values of both 
Ua and have been taken as 0.1. This is smaller than the maximal values of 0.87 and 0.37 specific 
to the third-order upwind scheme for the advection and the second-order central discretization for 
diffusion, respectively, together with a two-stage Runge Kutta time integration method [32]. 

Then, starting on the finest grid level, the particle fluxes are computed on each grid and corrected 
using the flux correction formulas (4.8). Then the first stage of the Runge Kutta step (3.18) is carried 
out. The field induced by the density predictors ct" and is again computed using the procedure 
described in Sect. 4.4. The new boundary conditions for the particle densities are computed on all 
grids TY""'*^, and the second stage of the Runga-Kutta method is carried out in the same way as the 
first stage. We then obtain the density distributions 17"+^ and on the set of grids 'H"''^. 

Following the procedure described in Sect. 4.3, after restriction of fine grid values to the parent grids, 
the grids 7^"+^''^ for the next time step are computed using the refinement monitor (4.3). The densities 
are then mapped from the grids K"'*^ to the grids and the new boundary conditions are 

computed. We thus obtain the density distributions cr'*+^ and p""*"^ on the set of grids 7^"'+^''= at the 
new time r""*"^. 



4-6 Relations with other refinement algorithms 

As mentioned before, the initial attempts to solve the system (3.1)-(3.4) were done by using the 
existing adaptive finite difference code VLUGR [34]. This code failed for our problem. Another un- 
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successful attempt was made using the adaptive finite element code KARDOS [42]. Both these codes 
use implicit time stepping, and they do not take into account the unstable behavior of the solution 
ahead of the front. 

Since our algorithm uses explicit (Runge-Kutta) time stepping, the refinement procedure is relatively 
simple. In particular, the computation for the electric potential becomes decoupled from the density 
updates in the continuity equations. This allows for a tailored approach to these separate problems. 
The potential updates are performed using nested fast Poisson solvers [39], which requires rectangular 
(nested) regions for these sub-problems together with a high-order intcrpolant for a global error 
indicator. For the continuity equations, on the other hand, the computational regions can be chosen 
quite small, essentially confined to the streamer itself. For these equation the simple error monitor 
based on local curvature performs well. 

Local time stepping, for the different grids on parts of the domain, have not been considered in this 
paper. Although savings might be expected for the continuity equations, where we now use the same 
time step dictated by the finest mesh, such approach would lead to the complication that also the 
potential needs to be updated locally, and that charge conservation is no longer straightforward. 

A grid refinement approach with fixed, a priori chosen, grids but with local time stepping has been 
discussed in [43] for a, somewhat related, system of plasma equations. The main difference with our 
problem is that this system is considered with low electric fields but on a much larger time scale, 
with movement of ions described by Euler equations. The time step restriction (3.23) for dielectric 
relaxation then becomes very severe. To overcome this restriction, an implicit treatment of the po- 
tential in the electron drift fluxes is required. In the approach of [44,43] this is done through (3.16) 
together with a low-order prediction update for the electron densities. The other processes, including 
higher-order corrections, are solved in a time-split fashion. With local time steps, synchronization at 
the new (global) time level is needed to ensure the relevant conservation properties to hold, such as 
charge conservation. Apparently, the instabilities in the leading edge are not that much of a problem 
on the time scales considered in these papers, probably because the charged fronts are smoother. 

Some results based on non-uniform, moving grids have been reported [24], for the simulations of 
streamers originating from point electrodes. The regridding approach used in [24] allows only for a 
fixed number of grid points, and does not seem to take into account properly the leading edge. In 
contrast, the method presented in this article enables a a fine grid to be put wherever needed, in 
particular, over the leading edge. 

We are not aware of other grid refinement approaches for systems of plasma equations that do take 
this leading edge instability into account. Related problems have been reported however in [45] with 
moving mesh methods for the ID Fisher equation ut = Uxx + 7^(1 — u). To overcome instability 
of the numerical moving mesh scheme, a special monitor function was advised in [46]. This has 
essentially the effect of mesh refinement ahead of the front. It seems that extensions of this moving 
mesh approach to multidimensional reaction-diffusion equations or more complicated systems have 
not been implemented yet. 
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5 Tests of the algorithm on a planar ionization front 



The implementation of the grid refinements is first tested on the evolution of a planar streamer front. 
As analytical results for this case are available [12], any errors in the implementation can be identified. 
Moreover, conclusions on the tolerances and interpolations can be drawn for the more general two- 
dimensional simulations. 

A one-dimensional simulation is carried out with the two-dimensional code by letting initial and 
boundary conditions depend on z only and not on the radial coordinate r. Specifically, the initial 
ionization seed 

(T(r,z,0) = /9(r,z,0) = 10-2 e-^^-^i")" (5.1) 

is located at z = and a constant background electric field = — |£;,|ez is applied. The spatial 
region in this specific example is (r, z) G [0, L] x [0, L] with L = 1024. The boundary conditions for 
the electron density and the electric potential in this specific case are 

da , , do , , da , 

a(r,0,r)=0, — (r,L,r)=0, — (0,z,r)=0, — (L,z,r)=0, 

^ (5.2) 
0(r,O,r)=O, |^(r,L,r) = |56|, |^(0,^,r)=0, |^(L, 2;, r) = . 



In this situation, the electric field does not depend on r. It can be written as £(r,z,T) = £{z,T)e2,, 
and it can be obtained directly from the charge densities by integrating V ■ £ = p — a from Eq. (2.10) 
along the z-direction and using the boundary condition £ = — \£b\ for the electric field at z = L. The 
result is ^ 

£{z,r) = -\£b\+ [ {aiz',T)-p{z',T)) dz'. (5.3) 

J z 

This means that it is not necessary to calculate the electric potential ^. Rather the electron and 

positive ion density at time determine the electric field at each cell vertex by discretizing 
Eq. (5.3). Starting from the value at z = L, which corresponds to j = M on a grid with M grid points 
in the z-direction, we thus obtain: 

^M+i = -\£l>\ , = + Az(a- - p]) for j<M. (5.4) 

The electric field strength in the cell centers is then taken as the average field of the corresponding 
vertices. 

We emphasize that in this section, in which the particular case of a one-dimensional streamer front is 
considered, we use Eqs. (5.4)-(5.5) to compute the electric field, rather than the Poisson equation. This 

speeds up considerably the computations, and enables us to test the refinements for the continuity 
equations only. For tests on the refinements of the Poisson equation only we refer to [39]. 

Fig. 5 shows the temporal evolution of the initial ionization seed (5.1) located at z^ = 31, with a 
background field \£h\ = 1 and a diffusion coefficient D = 0.1. The results shown are obtained on a 
uniform grid with mesh spacings Ar = 32 and Az = 1/4, which will be considered afterwards as a 
reference solution. Zb is chosen large enough that boundary effects at 2; = do not matter, and that 
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Fig. 5. Numerical solutions of the planar streamer front with \E^\ = 1 and D = 0.1 for < r < 200 with 
time steps of 10. The solution is computed on a uniform grid with a mesh spacing of 1/4 in the direction 
of propagation. The thick lines correspond to multiples of 50 for r. 

the initial maximum is well presented even in the case treated later in which a coarse grid with mesh 
size Az = 2 is used. After an initial growth until time r ~ 25, an ionization front emerges that moves 
with about constant velocity to the right in the direction opposite to the background field. 



The front velocity as a function of time is plotted in Fig. 6. It is derived from the numerical dis- 
placement of the level = 10~^ within a time interval of 10. The front decelerates and eventually 
approaches a value somewhat below the asymptotic front velocity [12] 

V* = \£b\ + 2^ D\£b\e-'^/\^b\ (5.6) 

which for these particular values of the background electric field and the diffusion coefficient is equal 
to 1.3836. The numerical velocity at large times is around 1.365 which corresponds within an error 
margin of 1.5% to the asymtpotic value. We expect to obtain even better results on finer grids, but we 
focus in what follows on the performance of the refinement algorithm compared to this uniform grid 
computation. For further discussion of the results we refer to Appendix B. Moreover, deviations from 
this asymptotic value can be derived analytically for different numerical schemes, and that this is 
illustrated in Section 5.6.6 of [25] for an explicit and a semi-implicit time discretization of a diffusion 
equation. 



In the following illustrations we have computed the temporal evolution of the densities and the electric 
field on a fine {Az = 1/4) and on a coarse (Az = 2) Tiniform grid as well as on locally refined grids. In 
the latter case, the coarsest grid has a mesh spacing of Az^ = 2 which we refine up to a finest mesh 
width of Az-^ = 1/4, thereby allowing for three levels of refinement. The electric field is again computed 
using Eq. 5.4 rather than through the Poisson equation for the electric potential, which speeds up the 
computations. The refinement algorithm for the continuity equations is as explained in Sect. 4.3. To 
demonstrate that the leading edge has to be included in the refinement, results with the "standard" 
refinement (i.e. without including the leading edge) will be compared to those that do include the 
leading edge. In this leading edge of the ionization front the densities decay exponentially, but the 
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Fig. 6. Numerical front velocity (thin solid line) as a function of time, compared to the theoretical front 
velocity (thick line, derived in Appendix B) and the theoretical asymptotic front velocity v* = 1.3863 

electric field strength is such that the reaction term is non-negligible. From theoretical studies [12] as 
briefly recalled in Appendix B it is known that this region is very important since it determines the 
asymptotic dynamics of the front. 

For the present problem we use the a priori knowledge that the front moves to the right and the 
leading edge is the region ahead of the front. The standard criterion reads, 



refine all grid cells j where 



u = a, a — p 



(5.7) 



Second, we use the same criterion but now with the inclusion of the so-called leading edge in the 
refined region, taking into account the cut-off of densities below the continuum threshold of 10~^^, 
i.e.: 

refine all grid cells j obeying criterion (5.7), 

extend the refined region in the propagation direction (5-8) 
to all {zj\aj > 10-^2}. 



The upper and lower plot of Fig. 7 show the results of the one-dimensional streamer simulation com- 
puted with the criterion 5.7 and 5.8, respectively, together with a coarse and a fine uniform grid 
computation. The time step is such that the Courant numbers for both advection and diffusion are 
at most 0.1, which is sufficiently small to render the temporal errors negligible. 

The value of Zb in the initial pulse (5.1) has been chosen such that the maximum of this Gaussian 
pulse was situated in a coarse grid cell center. Moreover, if a grid refinement is needed at the very first 
time step, the initial condition on the finer grids is not interpolated from the coarser grid values, but 
is calculated directly from Eq. (5.1), so that the initial pulse on the finer grids is computed without 
interpolation errors. 

The results with a fine uniform grid can be considered as a reference solution, and it is clear that on 
the coarse uniform grid the front moves too fast and is too smooth. This is due to the large amount of 
numerical diffusion introduced by the coarse grid. As can be concluded directly from the expression 
for the asymptotic front velocity, a larger diffusion constant leads to a larger velocity, and a larger 
velocity makes the front smoother [12]. 

The same is observed on the grids refined according to the standard criterion, even with a very strict 
tolerance the front is badly captured. This also appeared with other standard refinements codes, e.g. 
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Fig. 7. Performance of the refinements without and with inclusion of the leading edge. Plotted is a{z,T) 
as a function of 2; for r = 200. Upper row. results on grids refined with the standard criterion with a 
tolerance of 10~^ (circles) and 10^^ (pluses) compared to a fine grid computation {Az = thick solid 
line) and a coarse grid computation (Az = 2, thick dash dotted line). Lower row: the same, but now with 
leading edge included in the refinement, only shown with a tolerance of 10^^. 



VLUGR [34] The major conclusion from this failure is that the refinement takes place at the wrong 
place, namely at the regions with steep gradients. This is not where the dynamics of a pulled front 
is determined. This occurs rather in the leading edge where any perturbation of the linearly unstable 
state will grow [12]. Therefore the method in which the leading edge is included in the refined region 
gives even with a relatively low tolerance much better results than the standard refinement strategy. 

In Table 1 the front velocities Vf are listed for the fine {Az = and the coarse {Az = 2) uniform grids, 
as well as for refined grids {Az'^ = 2, Az-I" = 3) with or without inclusion of the leading edge, and with 
linear or quadratic interpolation of the boundary conditions; the tolerance is set to = ^a-p = 10~^. 
The velocities in the table have all been determined from the displacement of the maximal electron 
density between r=250 and 262.5. As explained before [12], the numerically computed front velocities 
should be smaller than the theoretical value, which is indeed the case for the fine grids. 

We also looked at the results obtained by discretizing the advective term with a first-order upwind 
scheme. We concluded from those tests that this scheme performs very badly, quite in contrast to 
what is said in [24] . The amount of numerical diffusion introduced by that scheme completely changes 
the asymptotic velocities on realistic grids. Moreover, numerical diffusion can be expected to over- 
stabilize the numerical scheme [32] and to suppress thereby interesting features of the solutions such 
as streamer branching. 

From these tests it appears that the grid refinement based on a simple curvature monitor works well 
provided the leading edge is included in the regions of refinement. This is in accordance with [12] 
where the importance of the leading edge for dynamics of a planar front is discussed. 
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Table 1 

Front velocities Vf on various grids. The refinements have been carried out with = ^a-p = 10~^. The 
exact asymptotic velocity is v* = 1.3836. For comparison, results for the Ist-order upwind scheme on the 
fine grid are added. 



method 


Vf 


uniform A2; = ^, Ist-order upw. 


1.585 


uniform Az = |, flux limited 


1.365 


uniform Az = 2, flux limited 


1.448 


standard refinement 


1.469 


refinement including leading edge 


1.365 



6 Performance of the code on streamer propagation with axial symmetry 



A planar ionization front as treated in the previous section is mainly of theoretical interest. Genuine 
ionization fronts are curved around the streamer head, which leads to field enhancement ahead of 
the front. We now consider the streamer propagation with cylindrical symmetry (2D case), which 
differs substantially from the planar front (ID case), mainly due to the field enhancement ahead of 
a curved front. In particular, the electric field cannot be calculated as easily as in Eq. (5.3), but has 
to be computed through the Poisson equation for the electric potential. The tolerance e^-p will play 
a non-negligible role in that case. Also, as the electric field is enhanced immediately ahead of the 
ionization front and decreases further ahead, the leading edge region of maximal linear instability of 
the non-ionized state is bounded. 

In this section the refinement algorithm as presented in Sect. 4 will be applied to a streamer initiated 
by a Gaussian ionization seed situated on the axis of symmetry at the cathode, z = Q. 



a(r,z,0) =p(r,z,0) = lO'^exp {-L^^^ . (6.1) 

The computational domain is Lr=1024: and L^=2048, and the background electric field is set to 
\Sh\ = 0.5. We use homogeneous Neumann boundary conditions at the cathode (z = 0), which, 
together with the ionization seed placed on the cathode, means that electrons fiow into the system. 
For other boundary conditions we refer to Sect. 2.3. 

The temporal evolution of the ionization seed (6.1) under these conditions is shown in Fig. 8. A 
uniform grid with grid size Ar = Az = 1 covers the domain where both the electron and positive 
ion densities are strictly positive. This is a relatively coarse grid, but, as mentioned in Sect. 4.1, the 
use of finer grids would require too many grid points for the FiSHPAK routine to handle within an 
acceptable accuracy. Therefore we will use this uniform grid computation as a reference to test the 
performance of the adaptive refinement method. Fig. 8 shows the snapshots of the electron and net 
charge density distributions as well as the electric field, at r = 75, 150, 225 and 300. 

At T = 75, the maximal deviation of the electric field strength from its background value is around 
0.4%, and space charges do not play a significant role yet. This is the electron avalanche phase, during 
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Fig. 8. Propagation of a Gaussian ionization seed as given by Eq.(6.1) in a background field \E\)\ = 0.5. 
The snapshots correspond to times t=75, 150, 225 and 300. The upper row gives the logarithm of the 
electron density a, the middle row the net charge density a — p and the lower figure the electric field 
strength \E\. A uniform grid Ar = Az = 1 has been used for both the continuity equations and the 
Poisson equation. 
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which the Gaussian electron density distribution advects with an almost constant velocity, undergoes 
a diffusive widening, and grows due to impact ionization, leaving behind a small trail with positive 
ions [47]. 

At r = 150 the space charges concentrate in a layer around the streamer head, as can be seen in 
the second column of Fig. 8. Due to its curvature, this space charge layer focuses the electric field 
towards the axis of symmetry, thereby enhancing it. The body of the streamer is sufficiently ionized 
(and its conductivity therefore enhanced) and the electric field in the streamer body is lower than the 
background field. The space charge layer becomes thinner and denser in time, and the electric field is 
increasingly enhanced, as can be seen in the third column of Fig. 8. 

Eventually, the streamer becomes unstable, and branches. The beginning of this branching state is 
shown in the rightmost column of Fig. 8. 

Let us now consider the effect of the refinements on the streamer dynamics, as well as the effect of 
cutting off the densities that are below the continuum threshold. To this end, we run the simulations 
with the same parameters - background field, initial and boundary conditions - as in the previous 
subsection, but we now allow one level of refinement for the continuity equations, and seven levels of 
refinement for the Poisson equation. The finest grids in both cases have a finest mesh size Ar-'' = t^z^ = 
1. In the next section, more levels of refinement will be used. In Sect. 5 we could see that, provided 
the leading edge of the streamer front is included in the refinement, and a quadratic interpolation is 
used to provide the boundary conditions for finer grids, a tolerance for the continuity equations of 
Ea = 10~^ is well suited. Moreover, the error in the spatial discretization of the net charge density 
cr — p is induced by the discretization of the drift and diffusion terms in Eq.(2.8). Hence it is natural 
to take the tolerance for the net charge density equal to that for the electron density. 

The choice for the tolerance for the refinement of the Poisson equation is less straightforward. In one 
dimension, the error in the second order discretization of the Poisson equation (2.10) reads 



^'^ = ^8^ = ^^^- ^'-'^ 

Therefore, in the one-dimensional case, the curvature monitor for the net charge density will also give 
the error in the spatial discretization of cj). In higher dimensions however, this correspondence does 
not strictly hold anymore. Nevertheless we assume that the tolerance for a — p will still give a good 
estimation for the error in the solution of the Poisson equation, and we therefore take 



€0 = e^-p = = 10 ^ . (6.3) 

The number of grid points in the r— direction contained in each fine grid for the continuity equations 
has be chosen as Mq = 32. The net charge density distribution at r = 75, 150, 225 and 300 are plotted 
in the upper row of Fig. 9, together with the grids on which the continuity equations have been solved. 
The cquipotential lines are shown together with the grid distribution for the Poisson equation in the 
lower panel of Fig. 9. 

Since the streamer front is not planar anymore, it is necessary to investigate the effect of the refine- 
ments not only on the axis, but in the whole streamer. The evolution of the half maximum contours 
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of the net eharge density is depicted in the left plot of Fig. 10, both in the uniform grid computation 
without cut-off below the continuum threshold, and with the refinements including the cut-off. There 
is an excellent correspondence between the uniform grid computation and the one where all grids 
are refined. At some places there is a slight difference between the results, but this is also a conse- 
quence of the results being plotted for the coarsest grid, i.e. Ar*^ = Az^ = 1 for uniform grid, and on 
Ar'^ = Az^ = 2 for the refined grids. The grid points on which the results are plotted therefore do not 
coincide, and the contours consequently might be slightly different. Up to branching however, there 
is no significant error in either the radius or the width of the space charge layer. A minor effect of the 
refinement only becomes visible once the streamer has become unstable. A more detailed investigation 
of this branching will follow in a later paper. 

The leftmost plot of Fig. 10 has to be completed with the evolution of the axial charge density 

distribution of the streamer, in order to control not only the position of the half maximum contours, 
but also the maxima themselves. This is shown in the middle and rightmost plots of Fig. 10, where the 
axial charge density and electric field strength distributions, respectively, are shown for the same times 
as the leftmost plot of Fig. 10. Again it appears that the adaptive grid refinement method produces 
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Fig. 9. Grid distributions for the continuity equations (upper row), and for the Poisson equation (lower 
row). The times are the same as in Fig. 8. For the continuity grids, the black regions correspond to the 
coarsest grid with Ar'^ = Az^ = 2, and which covers the domain on which the particle densities are above 
the continuum threshold. The gray regions are covered with the fine grids with Ar = Az = 1. The two 
coarsest grids - with Ar = Az = 128 and 64 - on which the Poisson equation has been solved cover the 
whole computational domain, which is filled in black. The finer grids with cell size 32, 16, 8, 4 and 2 are 
plotted in lighter gray shades, the white grid corresponds to Ar^ = Az^ = 1. 
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Fig. 10. Comparison of uniform grid results with those obtained using the grid refinement strategy. Left 
panel: evolution of the half maximum contours of the net charge density. Middle panel: evolution of the 
net charge density on the axis of symmetry, (cr — p){0,z,t). Right panel: evolution of the electric field 
strength on the axis of symmetry, \£\(0, z,t). The results in all panels are shown for r = 25 to 325, with 
steps of 25. Solid line: uniform grid (Ar = Az = 1), dots: grid refinement (Ar = Az = 1, 2). 

results that coincide very well with the uniform grid computation. 

Finally, because of the charge conservation during an ionization event, it should be verified that 
the refinements are indeed mass conserving, as has been taken care of in Sect. 4.3. This is shown in 
Fig. 11, where we have compared a uniform grid computation where the densities below the continuum 
threshold have not been cut off, with the refined grid computation. It is clear that neither the cut-off, 
nor the refinement disturb the total number of particles in a visible way. 

The gain in memory obtained with the refinement method can be observed through the number of 
grid points that are used to solve the continuity equations and the Poisson equation. As can be seen 
in Fig. 12, the number of grid points used are one to three orders of magnitude smaller than in the 
computations as performed in [14], where a uniform grid covers the whole computational domain. For 
the gain in CPU time, we refer to table 2 in the next section. 

We emphasize that in the previous test the choice for two levels of refinements has been made in 
order to compare the results with uniform grid computations. In later simulations we will use more 
refinement levels, thereby reaching much smaller mesh sizes. We can however extrapolate the outcome 
of these tests to cases using more levels. 
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Fig. 11. Evolution of the total number of electrons (solid line) and net charge (dashed line). The lines 
correspond to a uniform grid computation without cutting off the densities below the continuum threshold, 
the dots to the refined grid computation. 
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Fig. 12. Number of grid points as a function of time. The thin dashed line represents the number of grid 
points corresponding to a grid of grid size Ar = Az = 1 over the whole domain, i.e. 1024x2048. The 
thick lines give the temporal evolution of the number of grid points when the refinement algorithm is 
applied up to a finest mesh size of Ar = Az = 1. Solid line: continuity equations, 1 level of refinement. 
Dashed line: Poisson equation, 6 levels of refinement. 

7 Accuracy requirements for the streamer simulations 

To illustrate the use of the above algorithm, we present results in a new parameter regime, namely 

in a very long gap with relatively low background field, and a short gap with high field. Wc present 
results on different finest mesh sizes, and investigate the convergence of the solution on decreasing 
the finest mesh. 



7. 1 Long streamers in a low electric field 

We consider a gas gap on which a background electric field \€b\ = 0.15 is applied. The negative 
electrode (cathode) is placed at z = 0, the positive one (the anode) at a distance = 2^^ =65 536. The 
radial boundary is situated at Lr = 2^15) =32 768. For N2 at atmospheric pressure, this corresponds 
approximately to an inter-electrode distance of 15 cm and an electric field of 30 kV/cm. 

The initial seed (2.13) is placed on the cathode at z = (25 = 0). The maximal initial density is 
ao = 1/4.7 and the e- folding radius of the seed is Ri, = 10, which correspond to a density of 14 cm~^ 
and a radius of 23 /um, respectively, for N2 under normal conditions. This relatively dense seed enables 
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us to bypass the avalanche phase. If we would start with a single electron at this value of the electric 
field, analytical results [47] show that the transition to streamer would occur after the avalanche has 
traveled a distance d « 18Q;~^(|£b|) « 14 000. Using a dense seed accelerates the transition time 
considerably. The dense seed mimics streamer emergence from a needle inserted into the cathode or 
from a laser induced ionization seed. 

The coarsest grid for the continuity equations has a mesh size of 64 in both directions, the coarsest 
one for the Poisson equation a size of 8192. We first present the results when the grids are refined up 
to a mesh size Ar^ = Az^ = 2. 

The upper panel of Fig. 13 shows the electron density distribution on a logarithmic scale. There are two 
regions in the streamer: a rather wide one with low electron densities, and one which is much narrower 
and very dense. The lower panel of Fig. 13 clearly shows the negatively charged layer and its effect 
on the equipotential lines. These are close to each other ahead of the streamer tip, which indicates 
an enhancement of the electric field. In the interior field the distance between the equipotential lines 
is slightly larger than outside the streamer, which implies that the field in the streamer body is 
somewhat smaller than the background field. 

Let us now compare these results with those obtained on coarser grids. We have run the simulations 
on grids that were refined up to mesh sizes of Ar^ = Az^ = 4 and Ar^ = Az^ = 8, thus the finest 
grids in these cases are twice respectively four times coarser than those on which the results shown 
in Fig. 13 have been obtained. 

Fig. 14 shows the influence of the mesh size by means of the net charge density distribution and the 
electric field. The leftmost plot of this figure shows the evolution of the half maximum contours of the 
net charge density. The evolution of the density distribution and the electric field strength are shown 
in the middle and rightmost plot, respectively. Up to r « 5000 we see that the coarse grid simulations 
{Az^ = 8) already give convergent results. After this time, the front starts to move somewhat too fast 
on this coarse grid. This is due to numerical diffusion. The simulation with Az^ = 4 gives the same 
results as those on a grid with size Az^ = 2 up to r ~ 7500, after which numerical diffusion again 
makes the front move somewhat too fast. We emphasize However, the influence of the grid becomes 
significant only around r 10000, when the streamer head becomes unstable. The effect of the grid 
size on the streamer instability will be considered in more detail in a future paper. 

We notice that the solution on a grid with mesh spacing 8 eventually lags behind the flne grid 
computation. This non-monotonous behavior is due to a competition between two opposite effects: 
the numerical diffusion tends to accelerate the field, but also to reduce the maximum of the charge 
distribution, thereby reducing the maximal electric field strength, and thus the propagation velocity 
of the streamer. 

In Fig. 15 we compare the evolution of the maximal net charge density on the axis of symmetry for 
different choices of the finest mesh size. It shows that, the coarser the grid, the more underestimated 
the maximal net charge density becomes. This indicates that it is indeed numerical diffusion which 
smears out the space charge layer, thereby accelerating it and decreasing its density. 

We conclude that the simulations run on a finest grid with mesh size 8 will give non-negligible errors 
in the results. However, one more level of refinement, leading to a finest grid with mesh size 4, 
already gives acceptable solutions during the streamer regime. Moreover, from these results, we can 



33 



T=2500 t;=5000 t=7500 t= 10000 
5000| ■ 1 I ■ 1 I ■ 1 I ■ 




1000 2000 1000 2000 1000 2000 1000 2000 
r r r r 




1000 2000 1000 2000 1000 2000 1000 2000 
r r r r 



Fig. 13. Evolution of a streamer and the computational grids in a low background electric field \£h\ = 0.15. 
Upper panel: logarithm of the electron density with the grids for the continuity equations. The coarsest 
grid with mesh size 64, covering the domain on which the continuum approximation for the densities holds, 
is filled in black. The finest grids with mesh size 2 are white. The particle densities in the white region 
ahead of the streamer are below the continuum threshold and therefore not solved. Lower panel: net charge 
density (thick lines) and equipotential lines (thin lines) together with the grids for the Poisson equation. 
The white domains are covered with the finest mesh size of 2, the black regions are covered by a grid with 
mesh size 64. The gray region ahead of the streamer at t=2500 is covered by a grid with mesh size 128. 
The coarser grids (with mesh sizes up to 8192) are not shown. 
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Fig. 14. Influence on the mesh size on the numerical solution for a background field of 0.15. The left panel 
shows the evolution of the half maximum contours of the net charge density a — p. To show the structure 
clearly, the aspect ratio is not equal. The middle and the right panels show the evolution of the net charge 
density and the electric field strength on the axis, respectively. The times shown go from 500 to 11000 
with equidistant time steps of 500. The thick lines correspond to the times shown in Fig. 13. Solid line: 
Az-I'=2; dashed line: Az^ =A; dotted line: Azf=8. 




10000 



12000 



Fig. 15. Temporal evolution of the maximal net charge density on the axis. Solid line: Az^—l; dashed line: 
Az-I'=A; dotted line: Az^=8. 
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extrapolate that refining the grids even more (e.g. up to a finest mesh size of 1) will not lead to a 
significant correction of the results, and that the computational work required to run the simulations 
on such fine grids is not in proportion with the gain in accuracy. 



1.2 Fine grid computations at a high electric field 

We now consider the evolution of a streamer in a short overvolted gap. The inter-electrode distance is 
set to 2048, corresponding to 5 mm for N2 at atmospheric pressure. The background electric field is 
set to £h = — 0.5e^, which corresponds to 100 kV/cm. The initial seed (2.13) is placed on the cathode 
(z(, = 0). The maximal density of the initial seed is set to ctq = 10~^, and its radius to Ri, = 10. 

The evolution of the streamer with previous initial and boundary conditions, during the non-linear 
phase is shown in Fig. 16. These results have been obtained on finest grids with a mesh size of 
Az-I" = 1/8 for both the continuity and the Poisson equations. The coarsest grid for the continuity 
equation has a mesh size of Az'^ = 2, the one for the Poisson equation has a size of Az'^ = 128. 

The discharge clearly exhibits the streamer features which are, as in the low field case, a thin negatively 
charged layer that enhances the electric field ahead of it, and partially screens the interior of the 
streamer body from the background electric field. 

In order to investigate the convergence of the solution under decreasing grid size, we have also run 
the simulations on finest mesh sizes Az^ = 1/4, Az-^ = 1/2 and Az^ = 1. Fig. 17 shows the evolution 
of the half maximum contours of the net charge density, as well as that of the axial distribution of 
the net charge density and the electric field strength. 

We see that up to r = 200, a grid size Az^ = 1 gives the same results as a grid size Az^ = 1/8. 
After that time, the predicted front velocity on a grid size of 1 is much faster than that on the finer 
grids, whereas the maximal electric field does not differ that strongly from fine grid computations. 
This implies that the front propagates faster because of numerical diffusion. 

The finer grids on the other hand all give similar solutions up to r 300. After this time, the 
streamer head has become unstable, and the numerical details will influence the branching behavior 
of the streamer. Indeed, for Az-^ = 1/8, the streamer exhibits two branches that propagate off-axis, 
whereas for Az^ = 1/4 and 1/2 one branch continues propagating along the axis while the other does 
not. Therefore the results at the time of branching diflFer. An off-axis branching results in a decrease 
of the maximal fleld and densities on the axis, whereas an on-axis branching on the contrary results 
in an increase of these quantities. This is why the maximal net charge density and field strength on 
axis are smaller in the Az^ = 1/8 case than in the Az^ = 1/4 and 1/2 cases. 

We emphasize that the different branching form is due to the instable nature of the streamer: the 
streamer approaches a bifurcation point, and the further evolution is set by minor details, like the 
numerical grid, in a chaotic manner. There is, however, convergence of the branching time, which shows 
that unlike the behavior after branching, the onset of the instability is not triggered numerically. 

In Table 2 we summarize the efficiency and accuracy of the refinement algorithm in terms of the CPU 
time, memory usage and accuracy in the spatial electron distribution for different finest mesh spacings 
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Fig. 16. Evolution of a streamer and the computational grids in a high background electric field = 0.5. 
Upper panel: logarithm of the electron density together with the grids for the continuity equations. The 
coarsest grid has a mesh size Az'- = 2 and is filled in black. The finest grid has a mesh size of Az^ = 1/8 
and is white. The white region surrounding the coarsest grid has not been covered by a computational grid 
since the densities are below the continuum threshold there. Lower panel: net charge density (thick lines) 
and equipotential lines (thin lines) together with the grids for the Poisson equation. The black region is 
covered with a grid with mesh size 2, coarser grids are not shown on this scale. The grids are refined up 
to a mesh of Az-^ = 1/8. However, this finest grid is only used at t=375, in a very small part of the 
computational domain. The finest grid for the Poisson equation at t=100 has a mesh size of 1. At t—200 
and 300 the finest grid has a mesh size of 1/4. 
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Fig. 17. Influence on the mesh size on the numerical solution for a background field of 0.15. The left 
panel shows the evolution of the half maximum contours of the net charge density a — p. The middle 
and the right panels show the evolution of the net charge density and the electric field strength on the 
axis, respectively. The times shown go from 25 to 375 with equidistant time steps of 25. The thick lines 
correspond to the times shown in Fig. 16. Solid line: Az^—1/8; dashed line: Az-'^=l/4; dash-dotted line: 
Azf=l/2, dotted line: Azf=l. 



and refinement levels. To this end we start the simulations from a well developed streamer and run 
the different cases over a time T = 10, with a time step Ar = 5 • 10~^. This time step is small enough 
to get negligible temporal errors. The initial particle distribution is the numerical solution at r = 200 
obtained on a hierarchy of grids with a finest and coarsest mesh widths A2;'^=l/8 and Az^=2 for the 
continuity equations (see Fig. 16). For the Poisson equation, the coarsest mesh width was set to 64, 
and the finest on to 1/8. The solution at r = 210 with these grid settings is taken as the reference 
solution for the computation of the errors. This benchmark has been performed on the node of linux 
cluster, a 32 bit Opteron 1.4 Ghz with 16 Gb memory. 

The cases with refinements all use the same coarsest mesh spacings and refinement tolerance, the finest 
mesh spacing varying from Azf = 1/8 to Az-^ = 1. To illustrate the performance of the algorithm, 
we show results on uniform grids with spacing 1 and 2. A solution on a uniform grid of spacing 1/2 
could not be obtained because the Poisson solver would become very inaccurate (see Appendix A). 
We further note than on a regular pe (rather than the cluster node with the large amount of memory 
used for the benchmark) memory would also become a limiting factor. 

Since the grid configuration, and therefore the CPU time used per time step, is not constant, the 
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Table 2 

Performance of the refinement algorithm in terms of CPU time, memory usage and errors in the solution. 
Shown are results on refined grids (four upper rows), and on uniform grids (three lower rows). Az^ is 
the finest allowed mesh spacing for both the continuity and the Poisson equations, n^e^ is the number 
of levels for the continuity equation, niev = 1 corresponding to a uniform grid computation, for both 
Poisson and continuity equations. No- and are the number of grid points used for the continuity and 
the Poisson equations, respectively. The reference solution for the errors computation is the one obtained 
with Az^—\/^. An * denotes a case which could not be tested because of the inaccuracy of the Poisson 
solver and/or the lack of computational memory. 





niev 


CPU time (s) 




N4, 


l|e||Li 


l|e||L2 


l|e||L<^ 


1/8 


5 


11.06 


657 856 


93 584 








1/4 


4 


9.68 


193 024 


93 584 


5.09 • 10-9 


1.18 • 10-^ 


6.13 • 10-4 


1/2 


3 


6.67 


76 160 


69 452 


2.71 • 10"^ 


7.41 • 10-6 


4.62 • 10-3 


1 


2 


4.58 


21 632 


44 248 


1.04 • 10-^ 


2.78 • 10-^ 


1.64 • 10-2 


1/2 


1 


* 


* 


* 


* 


* 


* 


1 


1 


7.24 


2 097 152 


2 097 152 


1.54 • 10-^ 


2.78 • 10-^ 


1.26 • 10-2 


2 


1 


1.70 


524 288 


524 288 


6.19 • 10-^ 


1.77 • 10-4 


1.01 • 10-1 



average values over the run of the number of grid points and the CPU time are shown. The accuracy 
of the method can be characterized by the discrete Li, L2 and Loo-norms of the errors. For a grid 
function e = (e^ ) on a rrir x niz grid these norms are defined as 



||e||oo = max(|eij|). (71) 

2, after having restricted the values from 

Table 2 shows that, although the order of the error does not show a clear asymptotic behavior, there is 
an obvious decrease of the errors with finer mesh widths. We notice that the maximum of a at r=210 
on the coarse grid is approximately 2. The maximum relative error therefore becomes negligibly small 
for mesh spacings of 1/2 or smaller. For coarser grids the error however is significant, showing that 
the parameters used here do require a higher accuracy than that of 2 used in [13] or 1 used in [14]. 

Notice that the number of grid points A'^^ used to compute the electric potential is the same for 
Az-^ =1/4 and 1/8, which implies that in fact no grid with mesh spacing of 1/8 is used in the latter 
case. Apparently, a mesh width of 1/4 is sufficient for an accurate solution of the Poisson equation 
with this choice of the tolerance. 

The uniform grid calculations clearly illustrate the gain in efficiency of the refinements. Two orders 
of magnitude are gained for the computational memory. The gain of in computation time is less 
accentuated. This is due to the time gained by using less grid points being partially spent in the 
refinement procedure. On a regular pc with 1 Gb of memory, the gain in CPU is a factor 2.5 more 



'■,3 V 



The errors are computed on the coarsest grid level Az"^ 
finer grids where needed using Eq. 4.1. 
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pronounced. This is due to limitations in memory, due to which the computer has to swap, which 
slows down considerably the computations. Finally, the errors for Az-^ are of the same order in the 
uniform and in the refined case, which confirms the good performance of the refinement algorithm. 

A uniform grid with mesh spacing less than 1 could not be obtained for these specific simulation 
parameters. However, the errors with refinements show that finer grids are necessary for the solution 
to be reasonably accurate. We can therefore conclude that this refinement procedure allows us to gain 
computational time, but, more importantly, to gain so much computational memory that we are now 
able to reach a sufficient accuracy. 



8 SummEiry and conclusions 

We have presented an adaptive grid refinement strategy for the computation of negative streamers 
within the minimal model in a three dimensional geometry with cylindrical symmetry. The equations 
are rewritten in dimensionless quantities allowing the transcription of the results to arbitrary gases 
of arbitrary pressure. 

The numerical discretization arc based on finite volume methods. It uses a second order central scheme 
for the electron diffusion, and a flux limiting scheme for the advcction, so that the numerical diffusion 
is reduced and no spurious oscillations are introduced. The Poisson equation is discretized with a 
second order central scheme. The time stepping is achieved with an explicit two-stage Runge-Kutta 
method. 

The explicit time stepping method allows us to refine separately the grids for the continuity equations 
and those for the Poisson equation. The refinement criterion of the continuity equations is based on a 
curvature monitor of the solution, while that of the Poisson equation is based on an error estimation. 
It results in two series of nested grids, one for the continuity and one for the Poisson equation, that 
communicate with each other using adequate restrictions and prolongations of the densities and the 
electric field. The refinement method has been implemented in such a way that mass, and therefore 
charge, is conserved during the refinements. 

Tests on planar ionization fronts show that the leading edge has to be included in the refinements to 
capture the front velocity well. This is because the front penetrates a non-ionized high field region 
that is linearly unstable against even infinitesimal electron densities; in fact, the ionization front is 
a so-called pulled front whose velocity is determined in the linear leading edge region, and not in 
the nonlinear high gradient region of the front. A test on a genuine streamer emerging from a small 
Gaussian ionization seed in a relatively high background electric field shows that our refinement 
including the leading edge region performs very well. The space charge layer with its steep spatial 
gradients is also captured well by our refinements. Moreover the requested compTitational memory is 
three orders of magnitude smaller than in a uniform grid computation as performed in [14]. 

The algorithm enables us to investigate streamers with sufficient accuracy in new parameter regimes, 
namely in larger systems and /or in higher electric fields than previously. The results show how negative 
streamers increasingly enhance the field at their tip, both in lower and in higher background fields. 
The spatial gradients increase with the field [12] and therefore require an increasing spatial accuracy. 
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For simulating streamers in a background electric field as high as 0.5 in dimensionless units, we 
here used grids with meshes as fine as 1/8, much finer than the grids with meshes of 2 and 1 used 
previously [13,14]. Earlier simulations of negative streamers were carried out in a background field of 

0. 25, with a grid size of at least 2 [28]. These simulations show a field enhancement up to a value of 

1, and our tests show that a finer grid size is required for reliable results. 

In fact, since the characteristic scale of the inner structure of the streamer front is set by the ionization 
length a{\£\max) [12,17], one should take care that the grid size is fine enough to capture this length 
scale. As a rule, we suggest that the grid size should be at least four times smaller than the ionization 
length in the maximal electric field. For the case of a background field of 0.15, where the enhanced 
field grows up to 0.4, this implies a grid size not exceed 3 (where we used 2). For a background field 
of 0.5, the enhanced field grows up to 2; the rule then implies that the grid size should not exceed 
1/2 (where we actually used 1/4 and 1/8). 

Summarizing, we now have an efficient and reliable tool for simulating negative streamers within the 
minimal model up to the moment of branching which will be exploited in future studies of streamer 
physics. Additional effects like electron attachment or photoionization are important in more complex 
gases like air, in particular, for positive streamers. Such effects can be implemented in continuum 
approximation along the same lines. 

We finish with an outlook on open problems. First, fully three-dimensional simulations are required 
to follow streamer evolution after branching. We actually expect the branching time in cylindrical 
symmetry to be a close upper bound for the actual branching time in the fully three-dimensional 
situation [48,20]. Second, numerical solutions of our deterministic fluid model show that it is actually 
admissible for streamer propagation to neglect densities below the threshold of 1 particle per mm^ 
where the continuum approximation definitely ceases to hold; we have used this threshold on our 
computations with refinement. In this region which actually belongs to the leading edge of the pulled 
ionization front, the discrete nature of particles should be taken into account, such work is presently 
in progress. 



Appendix A: Testing the FiSHPAK fast Poisson solver 

The FiSHPAK routine is used to solve the Poisson equation. It is a fast Poisson solver but has limitations 
with regards to the number of grid points. We illustrate the instabilities with the FiSHPAK routine on 
large grids by an example. Consider the Laplace equation in a radially symmetric coordinate system 
{r,z) G {0,Lr) X (0,L,), 

' VV = ( - + AA^r^ + {z- iL,)2)) e-^{r'+{z-hL.r) , 

< (^(r,0) = 0(r,L^) = e-^('''+(^-5^-)'), (8.1) 
with Lr = Lz = 1, giving the analytical solution 4>{r, z) = e"'^^^^^^^'"^^/^^^^ 



41 




Fig. 19. The Li-errors (solid), L2-errors (dashed) and Loo-errors (dash-dotted) for \£\ in Eq.(8.1) on 
m X m grids, as a function of m. Upper panel, single precision; lower panel, double precision. 

The accuracy of the method can be characterized by the discrete Li, L2 and Loo-norms of the errors 
given in Eq. (7.1). Figure 18 shows tliese Lp-norms of the numerical error of the results obtained 
with the FiSHPAK routine on a (m x m)-grid, as a function of m, with A = 100. The upper and lower 
panel show the error of the single respectively double precision computations. In a similar way we 
determined the Lp-norms of the numerical error in the electric field strength \£\ij, computed with 
Eq. (3.15). They are shown in Figure 19. 

Up to a number of grid points of approximately 2000, the errors in the electric potential are of second 
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order, in agreement with the discretization. The errors in the field are also of second order, even 
though, at first sight, a first-order behavior might be expected, since it's the derivative of a quantity 
of second order accuracy. This can be explained through an asymptotic error expansion, which, for 
simplicity, will be carried in one-dimension. The second order discretization for will give 



= (t){zj) + /^z^v{zj) + O(Az^) , (8.2) 

with a principal error function v which will be smooth if the solution (j) is so. By considering local 
truncation errors it is seen that v should satisfy d^zV = —d*(p/12 with corresponding homogeneous 
boundary conditions. Then the discretized value of the field at the cell boundary becomes: 

i+l = = + - + 0{Az ), (8.3) 

which, using a Taylor expansion around , i, yields 



£j+^ = -OM^j+i) + 0{Az^) - Az^dzv{zj) + 0(Az3) = £{z.^i) + 0{Az^) (8.4) 

It is clear that the performance of the FiSHPAK routine decreases dramatically when more than ap- 
proximately 2000 grid points (in double precision) are used. We note that these numerical instabilities 
also show up on a x niz grid if either nir or niz are larger than 2000, approximately, in double 
precision. 



Appendix B: Analytical expressions for the planar front velocity 



The asymptotic velocity (5.6) and the convergence towards it can be derived analytically, following 
arguments in [25]. First, the analysis in [12] shows that a negative streamer front is a pulled front 
whose dynamics is determined in the leading edge. Linearizing the one- dimensional streamer equations 
in the leading edge around the field £b < ahead of the front gives 

where f{\£b\) = l^fel exp(— l/|£^b|). As the nonlinear front region has to be "pulled along" by the 
evolution of the linear perturbations, the evolution of the linearized equation (8.5) yields an upper 
bound for the velocity at each instant of time. 

For an initial condition of Gaussian form 

a{z,0) = (To exp (- \ , (8.6) 
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as used in the paper, the hnearized equation (8.5) is solved through 



.(..)=.oexp J^^exp(-fc^!^) (8.7) 



_ Rl + ADt ^""^ V Rl + "^Dt 

for all times r > 0. Solving this equation for z, the position Zf{T) of a fixed electron density level af 
is 



z±(r) = Zb + \Sb\T ±^Rl + ADt ^ f{\Sb\)T - 



T>1 



^ Zb + {\£b\ ± V^JiM) r + 0{VT). (8.10) 
For the pulled front moving to the right, we need z^{t)- The velocity of the level cr/ is 



v^Pir) = drzjir). (8.11) 

This equation determines the asymptotic velocity (5.6) that is actually independent of the parameters 
(Tq, Rb and zi, of the initial conditions. Moreover the time dependent velocity t;j"(r) is an upper bound 
for the actual velocity w/(t) of the full nonlinear problem. The two velocities are compared in Fig. 6. 
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